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Abstract 

Beamforming in ultrasound imaging has significant impact on the quality of the final image, controlling its 
resolution and contrast. Despite its low spatial resolution and contrast, delay-and-sum is still extensively used nowa¬ 
days in clinical applications, due to its real-time capabilities. The most common alternatives are minimum variance 
method and its variants, which overcome the drawbacks of delay-and-sum, at the cost of higher computational 
complexity that limits its utilization in real-time applications. 

In this paper, we propose to perform beamforming in ultrasound imaging through a regularized inverse problem 
based on a linear model relating the reflected echoes to the signal to be recovered. Our approach presents two major 
advantages: i) its flexibility in the choice of statistical assumptions on the signal to be beamformed (Laplacian and 
Gaussian statistics are tested herein) and ii) its robustness to a reduced number of pulse emissions. The proposed 
framework is flexible and allows for choosing the right trade-off between noise suppression and sharpness of the 
resulted image. We illustrate the performance of our approach on both simulated and experimental data, with in 
vivo examples of carotid and thyroid. Compared to delay-and-sum, minimimum variance and two other recently 
published beamforming techniques, our method offers better spatial resolution, respectively contrast, when using 
Laplacian and Gaussian priors. 


Index Terms 

Adaptive beamforming, linear inverse problems, beamspace processing. Basis Pursuit, Least Squares 


I. Introduction 


U ltrasound (US) imaging is one of the most fast-developing medical imaging techniques, allowing 
non-invasive and ultra-high frame rate procedures at reduced costs. Cardiac, abdominal, fetal, and 
breast imaging are some of the applications where it is extensively used as diagnostic tool. The new 
advances in beam formation, signal processing, and image display enlarge the US imaging potential to 
other fields like brain surgery, or skin imaging (e.g. |[T| and Q). 

In a classical US scanning process, short acoustic pulses are transmitted through the region-of-interest 
(ROI) of the human body. The backscattered echo signals, also called raw radiofrequency (RF) data, 
are then processed for creating RF beamformed lines. Beamforming (BF) plays a key role in US image 
formation, influencing the resolution and the contrast of final image. The most used BF method is the 
standard delay-and-sum (DAS) which consists in delaying and weighting the reflected echoes before 
averaging them. So far, its simplicity and real-time capabilities make it extensively used in ultrasound 
scanners. However, its weights are independent on the echo signals, resulting in beamformed signals 
with a wide mainlobe width and high side-lobe level. Consequently, the resolution and the contrast of 
final image are relatively low Q. Several adaptive beamformers (with weights dependent on data) from 
array processing literature were applied to US, the most common being the Capon or minimum variance 
(MV) BF 0. It offers a very good interference rejection and better resolution than DAS, allowing higher 
contrast ||5|. However, this method uses an estimated covariance matrix of the data and its main issue is 
the high computation complexity due to the calculation of the inverse covariance matrix. To overcome 
this, many improved versions of MV have been recently proposed (e.g. Q and Q), but still not adequate 
for real-time applications. In practice, in order to provide well-conditioned covariance matrices, diagonal 
loading, time and spatial averaging approaches were investigated, see e.g. 0. respectively Q. 

Recently, to improve the MV BF, Nilsen et. al. proposed a beamspace adaptive beamformer, BS-Capon, 


and unlike MV BF, they based their BF method on orthogonal beams formed in different directions |10|. 
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This technique was also applied by Jensen et al. to develop an adaptive beamformer based on multibeam 
covariance matrices GD. called multi-beam Capon beamformer. In their works, a covariance matrix is 
calculated for each range in the image, based on the idea that the beams were transmitted with different 
angles. Thus, the authors of [ fTT] | were able to reduce the computation time of MV BF, while improving 
the resolution of the point-like reflectors. Following a similar idea, Jensen and Austeng applied to US 
imaging a method proposed initially by Yardibi et al. 0. called the iterative adaptive approach (lAA). 
They obtained better defined cyst-like structures compared to conventional DAS, and better rendering than 
MV. 

The work presented here uses a similar idea of beamforming range by range. However, inspired from 
the source localization problems, we represent, for each range, the BF as a linear direct model relating the 
raw samples to the desired lateral profile of the RF image to be beamformed. This formalism allows us 
to invert the problem by imposing standard regularizations such as the £i or £ 2 -norms. These choices are 
motivated by the existing works in US image enhancement, that are typically based on Laplacian (e.g., 
p4| |) or Gaussian (e.g., | fT5] |) priors. Thus, the major contribution of this paper is the improvement of 
the existing beamforming techniques by combining the proposed direct model formulated in the lateral 
direction of the images with a regularized inversion approach. Moreover, we incorporated the proposed 
method with a beamspace processing technique, in order to highly reduce the number of the required US 
emissions. 

In contrast to existing BF methods in US imaging using regularized inverse problem approaches (see 
e.g., p^-[|^), our method does not use the system point spread function (PSF) in the direct model or in 
the inversion process. Thus, the proposed BF technique does not require any experimental measurement 
(e.g., [ |22] |) or estimation of the PSF (e.g., [ [2^ ) [ [T4| |, [|24]|). 

Laplacian and Gaussian statistics, two of the most common regularizations in such imaging problems 
(including US imaging applications such as deconvolution), are considered herein, allowing the reader 
to observe their influence on the results. Furthermore, our method opens new tracks for more complex 
regularization terms (e.g. [|25]|-[|T7]|) to further improve the results. The proposed approaches, generically 
named Basis Pursuit beamforming (BP BF), respectively Least Squares beamforming (LS BF) in this 
paper, were evaluated using different Field II simulated data and in vivo carotid and thyroid experimental 
data. Finally we compared our BF techniques with four existing beamformers: the conventional DAS, 
MV, multi-beam Capon, and lAA. 

The reminder of the paper is organized as follows. First, in Section we summarize the theoretical 
background of BF. In Section |n^ we describe the proposed BF method from a regularized inverse problem 
perspective. Details about the experiments and the results are given in Section |V| and Section [VT| concludes 
the paper. 


11. Background 

The main elements used to model the BF process are depicted in Fig. [T| We consider, without loss of 
generality, the particular setup of an M-element US probe (M can be also the number of active elements 
of the probe), with the transducer’s elements denoted by Um, with m = 0, • • • M — 1. We consider a trivial 
change of variable, such that the position of the m-th element is: 

Pm = - (M - l)/2](A/2), m = 0,1, • • • ,M - 1, (1) 

with the probe’s elements positioned symmetrically around origin. We have considered, as example, the 
pitch equal to A/2 (the spacing between elements is half of the wavelength A = c//o, c denoting the 
speed of the sound through soft tissue, and /o the transducer’s center frequency). 

A series of K focused beams are transmitted with different incident angles 9k, k = The 

returning echoes are recorded using the same US probe, being time-delayed, such that the time-of-flight 
is compensated, so the backscatter from the point of interest is summed up coherently. If we consider 
that each of the recorded raw signal after the time-delay compensation has N time samples, the size of 
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M-element US probe 



Fig. 1: The main US imaging elements used to adapt the array processing beamforming methods discussed 
in Section and Section III to medical US imaging. 


the recorded data from all the K directions will be M x N x K. Finally, the final RF US image is a 
collection of RF beamformed lines, each of which being the result of beamforming the raw RF signals 
coming from an emission in the direction 9k, k G {1,... ,K}, using M elements of the tranducer. 

The classical DAS BF can be expressed as: 

M 

Skin) = Amin)), n = l,---,N, 

m=l ^ ^ 

k = l,--- ,K, 

where Amin) is the time delay for focusing at the point of interest sample, being dependent on the 
distance between the m-th element and the point of interest, Wm are the BF weights, ym^ G is the 

raw data received by the m-th element of the US probe, corresponding to the emission steered at angle 
dk. A simplified form of (|^ can be formulated as: 


Sk = w^yk, (3) 

where G is the time-compensated version of ym'^ in @ for the k-th emission (for the sake of 

generality, we consider to be complex-valued data), w is the vector of the beamformer weights of size 
M X 1, and (•)^ represents the conjugate transpose. DAS BF selects the weights independent on data, 
solving: 

minw^w, such that w^l = 1, (4) 


where 1 is a length M column-vector of ones since the raw data was focused using time-delays. The 
solution of Q is: 

^ ( 5 ) 


If we replace Q in @ we get: 


wdas - 


h = ^1 yk, 


( 6 ) 


where l denotes the transpose. A common technique used in US is to apply weighting functions such 
as Hanning, or Hamming apodizations to @ to further reduce the sidelobes of s^, resulting in improved 
contrast of the beamformed image, at the cost of a slight lateral spatial degradation. 

Further details about adaptive BF in US imaging (i.e., the methods used for comparision purpose) and 
about beamspace processing are provided in Appendices respectively in Appendix 
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Fig. 2: The elements used to form the proposed model. 


III. Proposed method: Beameorming through regularized inverse problems 
A. Model formulation 

The main elements used to model the proposed method are depicted in the Fig. For sake of simplicity, 
let us focus our problem at a time-sample (range) n. The proposed BF method is sequentially applied in 
the same manner to each range. If G is the raw data after the compensation of the time-of- 


flight for the /c-th steering direction 9k, we can form the steering vectors as in ( |2^ . Let A be the M x K 
steering matrix containing the steering vectors in (|26l) for all 9k directions, k=l, 


,K: 


A. •) ^02 1***5 ^0h 


(7) 


Note that A is known and depends on the positions of the probe elements and on the locations to 
beamform. Thus, it is independent on the actual positions of the reflectors. 

For each range n, we want to estimate the signal corresponding to a reflector as a function of its location, 
that will contain dominant peaks at reflector positions. Thus, the main difference from the multi-beam 
Capon beamforming method is that instead of calculating the values of the weights as in ( |^ , that 
are further used to calculate the reflector’s signal, we are directly estimating the corresponding signal by 
considering the raw data yk[n] as observations. In other words, we want to obtain an estimate of the 
reflected echo x[n\ G through the observations Unfortunately one difficulty arises: since 

Vkin] is corresponding to only one emission, modeling our problem using raw data as observations to 
estimate the reflectors, requires high computational cost, since we are dealing with multiple directions. 
We recall that the size of raw data in our problem at a range n, is M x K. To overcome this issue, 
motivated by the results in [[28]| and [|29||, we propose to use the DAS beamformed data instead of the 


original raw data. In addition to data dimensionality reduction, it was shown in [ [28| | and p9| | that more 
accurate results may be achieved by proceeding in this way. Thus, we can formulate our model as follows: 


s[n] = A)x[n\ + g[n\, 


( 8 ) 


where s[n\ G is a lateral scanline of the DAS beamformed image formed as discussed in the Section 
|n| A is the steering matrix formed with Q, and g[n] an additive white Gaussian noise. If we denote by 
S the DAS beamformed image of size K x N, formed by juxtaposing the DAS RF lines Sk expressed in 
@ for all K directions, we consider s[n] the lateral scanline from S taken at the time-sample (range) n. 
Note that, since the transducers are emitting the same pulse, we assumed that x[n\ which is the unknown 


signal, is the same for all K emissions, and for all transmitters (see e.g. [30|). 













5 


The role of the multiplication of the steering matrix A with its conjugate transpose in ([8]) is to 
relate the position of the elements with the position of all K reflectors on a scanline. This relation is a 
result of considering on the one hand that the elements are impinging to the reflectors situated on the 
lateral scanline (the multiplication of A with a?[n]), while on the other hand the reflectors are impinging 
to the elements through their reflected pulses (the multiplication with A^). Hence, the result of the DAS 
beamformed scanline s[n] is related to the unknown signal x[n] through a direct linear model. Fig. 
offers a schematic representation of our model in ([8]). Thus, after the compensation of the time-of-flight, 
the received raw data, y[n\ at a range n is formed by multiplying the steering matrix A with the desired 
signal x[n], y[n] = Ax[n]. This multiplication could be sufficient for describing the proposed model if 
we are considering the raw data y[n\, as observations. However, since we are using the DAS beamformed 
data instead of the original raw data, we further take into account the geometrical relationship between 
the potential sources and the elements of the probe (through the multiplication with A^). 

B. Beamspace processing 

In order to solve we firstly apply beamspace processing, a common tool used in source localization 
approaches that reduces the computational complexity, while improving the resolution, and reducing the 
sensitivity to the position of the sensor (see e.g. p9] | and [|3T||). Its main purpose in US is to reduce the 
number of the US emissions, thus reducing the acquisition time and the computational load required by 
the BF process. We should note that our method of transforming the data into beamspace domain is totally 
different from the technique resumed in Appendix The main reason is that, by using the beamspace 
processing presented in p0| , we need all the acquired raw data for applying beamspace processing as 
described in ( |22| ). Hence, even if on one hand, the computational complexity required by the estimated 
covariance matrix inversion is reduced, on the other hand, it is increased by the operations required to 
transform the entire set of the raw data into its beamspaced correspondents. 

To overcome this, we based our idea on the beamspace processing techniques proposed in array 
processing (notably in source localization). More specifically, Malioutov et al. |j^ used a method that 
maps the data from full dimension space of the directions (DS) into a lower dimension beamspace (BS) 
through a linear transform prior to source localization processing. In our case, for each range n we project 
the data resulted by applying DAS BF, s[n], in BS before beamforming it through regularized inverse 
problems. To emphasize, s[n] E is projected on a lateral sampled grid of P « K locations. In other 
words, the proposed BF method, contrarily to all the other discussed BF methods, uses only P focused 
emitted beams among all the K transmissions to beamform a particular lateral scanline of K samples. 
Thus, the number of emissions is reduced by a factor of This will result in a reduced dimensionality 
of the data compared with the other BF methods, and an improved computational complexity compared 
with MV, multi-beam Capon, and lAA. 

Let z[n\ G be the beamspace transformed vector formed by sampling the DAS beamformed 

lateral scanline s[n] on a grid of P locations, see Fig. 

z[n] = (9) 

where D of size A x P is the beamspace decimation matrix, that will reduce the dimensionality of a 
vector from A x 1 to P x 1. Hence, since the decimation factor is D has all elements zero, except 
the elements dij with j = yi, that will get the value 1. Similarly, the beamspaced steering matrix A^g 
of size P X M is formed, composed of the P beamspaced steering vectors: 

= D^A^ (10) 

Concretely, we form Abs & ^mxp ^^king from A^ each ^-th steering vector. So, the model formed 
by @ after applying beamspace processing, becomes: 
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z[n] = D^s[n] = {ABsA)x[n] + g[n], (11) 

where x[n] of size K x 1 is the lateral profile at range n to be estimated. Thus, we can see o as 
an inverse problem, where z[n] is the DAS beamformed data corresponding io P < K emissions, and 
considered as the observation data. 


C. Beamforming through regularized inverse problems 

Given the ill-posedness of the direct model in ( [TT] ), we propose hereafter to invert it using standard 
regularization techniques. For achieving this, a cost function, denoted by J{x[n]), consisting into the 
linear combination of two terms is considered. The first term, denoted by Ji{x[n]), represents the data 
attachment, while the second, denoted by J 2 {x[n]), is the regularization prior: 


J{x[n]) = Ji{x[n]) + XJ 2 {x[n]), 


( 12 ) 


where A is a scalar, called regularization parameter, that adjusts the trade-off between the fidelity to the 
data and the regularization term. A large A will strongly favor the a priori about x[n\, while a small A 
gives a high confidence to the observations. Keeping in mind that the additional noise in ( fTTj ) is Gaussian, 
the data attachment term is expressed by an ^ 2 -norm, giving the following cost function: 

J{x[n]) = \ \z[n] - (AssA)x[n]\\l + XJ 2 {x[n]). (13) 


In this work, the choice of the regularization term J 2 {x[n]) is guided by the existing literature on 
statistical modelling of US images, previously applied to various applications such as image deconvolution 
or segmentation (e.g. [ fT4| |, [[32||). It has thus been shown that Laplacian and Gaussian statistics are well 
adapted to model US images. For this reason, we give in the two following paragraphs the mathematical 
derivations and beamforming results using ii norm (the sum of absolute difference) and £2 norm (or the 
Euclidean norm, that is the sum of squared difference) regularization terms. We should note that while 
the first will promote sparse solutions, the latter will promote smoother results. 

The choice of a quadratic data fidelity term is related to the additive zero-mean Gaussian assumption 
on the noise. We emphasize that the noise considered in our paper is different from the multiplicative 
speckle noise generally assumed to affect envelope images in ultrasound imaging. Instead, the additive 
noise considered in our paper affects the RF data and is caused by the acquisition process. The same 
model has been previously used by several authors (e.g. 0 and poll). 


1 ) Laplacian statistics through Basis Pursuit (BP): Considering that the signal x[n] to be beamformed 
follows Laplacian statistics, the minimization of the cost function in ( [T3] ) turns into the optimzation 
procedure in ( fid) ), usually called Basis Pursuit (BP) in the literature p^. 


XBp[n\ = argmin(|| 2 :[n] - {A^s^)^[n ]\\2 + 


(14) 


where || • ||i denotes the /i-norm. The mi nimization problem ( fT4l ) is convex, hence continuous, and has 
one global minima for any A > 0. 

Herein, we used the well known YALLl to solve ( fT4l ) p4| |, a software package that contains imple¬ 
mentation of alternating direction method (ADM), that solves also BP. A comparison of the six most used 


BP implementations is done in p5] | and three of them were also compared in p^ with application in 
underwater acoustics, where is shown that YALLl gives best performances for real-time applications. 
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TABLE I: Parameters of simulated and experimental images 


Parameters for simulation of: 

Point scatterers 

Fig.| 

Cyst simulation 

Fig.| 

Simulation of cardiac Image 

Fig.| 

Experimental carotid 

Fig.|^ 

In vivo thyroid 

Fig.[^ 

Transducer 

Transducer type 



Linear array 



Transducer element pitch [/urn] 

256 

256 

192.5 

no 

120 

Transducer element kerf [/xm] 

20 

20 

38.5 

25 


Transducer element height [mm] 

5 

5 

14 

4 


Central frequency, /o [MHz] 

3 

3 

4 

7 

12 

Sampling frequency, fs [MHz] 

100 

100 

40 

40 

40 

Speed of sound, c [m/s] 



1540 



Wavelength [/xm] 

513.3 

513.3 

385 

220 


Excitation pulse 



Two-cycle sinusoidal at /o 



Synthetic Aperture Emission 

Receive Apodization 



Hanning 



Number of transmitting elements 

64 

64 

64 

128 

128 

Number of receiving elements 

64 

64 

64 

128 

128 

Number of emissions (K) 

260 

260 

204 

192 

312 

The values of A for the simulated and experimental images 

BP 

0.5 

0.5 

0.2 

0.5 

5 

LS 

0.7 

1 

0.5 

1 

1 


2) Gaussian statistics through Least Squares (LS): To achieve smooth solutions of the proposed 
BE method, we modeled our problem with an £ 2 -norm based minimization function, and we used the 
Tikhonov regularized least-squared method (or rigid regression) for solving it [ |T7] |. The cost function is 
of the form: 


XLs[n] = aYgmm{\\z[n] - {ABsA)x[n]\\l + X\\x[n]\\l), (15) 

x[n] 

where || • ||2 denotes the f 2 -norm. For solving ( [TS] ) we used its analytical solution: 

XLs[n] = ({AssA)^{AssA) + XlKxK)~^iABsA)^z[n], (16) 

where Irxk denotes the identity matrix of size K x K. 

In order to obtain the BP and LS beamformed images, for each time sample n, with n e 
we estimate its corresponding lateral scanline XBp[n\ (using BP BE method) or (using LS BE 

method), and we are juxtaposing all the obtained scanlines, in the axial direction of the image. 

IV. Experiments 

For evaluating the proposed BP BE and LS BE approaches, we considered different types of simulated 
and experimental data. We compared our BE results with DAS (Section [II]), MV (Appendix |^, multi¬ 
beam Capon (Appendix |^, and lAA (Appendix BE methods. The simulations were made using the 
Field II program (see e.g. [ [38] | and ||^). The first simulation include a sparse medium, the second one 
contains a circular hypoechoic cyst in a medium with speckle, and the third one contains a simulation of 
the Short Axis (SAx) view of a cardiac image, as suggested in |j^. The first experimental data consists in 
a carotid that was recorded with an Ultrasonix MDP research platform. Finally, the second experimental 
data contains a thyroid medium with a malignant tumor. The thyroid data was recorded with the Sonoline 











Elegra clinical scanner, modified for research purpose. The parameters of the simulated and experimental 
data are presented in the Table |I| Note that for MV and multi-beam Capon beamformers, spatial and 
temporal averaging, as well as diagonal loading technique are used for the estimation of the covariance 
matrix, as discussed in the Section 

An important aspect is that, when applying the proposed BE methods, five times (^ = 5) less emissions 
were used in the beamforming process, by applying the beamspace processing presented in the Section 
III-B[ This hangs on for all the examples we are presenting in this paper. For these examples, reducing 


five times the US transmissions is optimal in terms of gain in resolution, while reducing computational 
time. 

The values of the regularization parameter A for all the presented examples are grouped in the Table |Ij 
The optimal value of A was chosen manually. We emphasize that this was the case for all hyperparam¬ 
eters of all comparative methods. Several studies exist in the literature for automatic estimation of the 
regularization hyperparameter (e.g. ED. p2| |, and [ |4^ ) that can improve the robustness of the proposed 
methods, at increased computational cost. Nonetheless, their implementation is beyond the scope of this 
paper. 


A. Parameters for the comparative methods 

The results of MV beamforming were obtained by using the implementation described in Q. The 
length of the spatial averaging window, L was defined as half of the number of the probe’s elements, i.e. 
L = ^. A temporal window of 10 samples was used in our examples and the diagonal loading parameter 
was fixed to A = The adaptive coherence method was applied to the MV BE method. 

The results of multibeam-Capon were obtained by using the multi-beam approach suggested in pT| |. The 
K emissions were uniformly distributed between ±30°. The beamspace transform down to 33 dimensions 
was applied, able to retain the variance for incoming narrowband far-field signals. The diagonal loading 
factor was set to 0.01. 


Finally, for lAA implementation we used the source code provided by the authors in pA] . The number 
of iterations was set to 15 for our examples. 

Note that for all the comparative methods, several parameters need to be carefully tuned in order to 
obtain acceptable results. However, using the proposed approach, only the regularization hyperparameter 
A needs to be set. 


B. Simulated point reflectors 

The medium contains 5 point reflectors, 4 of them aligned in pairs of 2 and separated by 4 mm, and 
the other laterally centered at 0 mm. They are located at axial depths ranging from 63 to 68 mm, with a 
transmit focus at 65 mm and a dynamic receive focalisation. 

C. Simulated phantom data 

To evaluate the accuracy, contrast and resolution of the aforementioned beamformers, a hypoechoic cyst 
of radius 5 mm, located at the depth 80 mm, in a speckle pattern. The speckle pattern contains 50000 
randomly placed scatterers, with Gaussian distributed amplitudes. This example was inspired from the 
simulation of a synthetic kidney example included in Field II software. The attenuation was not taken 
into account. 


D. Simulated cardiac image 

The SAx view is the cross-sectional view of the heart and is a well exploited perspective in echocar¬ 
diography, containing information about the left ventricle (LV) and right ventricle (RV). In our simulation 
we visualize the LV. The transmit focus point is set at 65 mm. The final image is ultra-realistic, the 
amplitudes being related to an in vivo cardiac image |j^. The number of scatterers was sufficiently large 
to produce fully developed speckle. 
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E. In vivo data: carotid 

The carotid ultrasound is a common procedure used to detect strokes or the risk of strokes due to the 
narrowing of the carotid arteries. The data was acquired from a healthy subject, with the Ultrasonix MDP 
research platform, attached with the parallel channel acquisition system, SonixDAQ. The linear ultrasonic 
probe L14-W/60 Prosonic® (Korea) of 128 elements was used. 

F. In vivo data: thyroid 

The thyroid ultrasound is done to visualize the thyroid gland to detect possible tumors or deformations. 
Two sets of data were acquired: first one, from a subject with a tumor and the other one from a healthy 
subject. For both acquisitions we used the clinical Sonoline Elegra ultrasound system modified for research 
purposes, and a 7.5L40 P/N 5260281-L0850 linear array transducer of Siemens Medical Systems, having 
the characteristics described in Table H 


G. Image quality measures 

Three image quality metrics were evaluated: the contrast-to-noise ratio (CNR), the signal-to-noise 
ratio (SNR), and the resolution gain (RG). They were computed based on the envelope-detected signals 
independent of image display range. 


Based on two regions Ri and R 2 belonging to two different structures, CNR is defined as [ [44] |: 

IhRi — hR2\ 


CNR = 


(17) 


cr 


Ri 


a 


R2 


where and are the mean values in the regions Ri, respectively R 2 , and and gr^ are the 
standard deviations of intensities in Ri, respectively R 2 . 

The SNR is defined as the ratio between the mean value /i and the standard deviation g in homogeneous 


regions [12|: 


SNR = 




(18) 


The RG is defined in [ [43] | as the ratio between the normalized autocorrelation function with values higher 
than 3 dB (computed for the DAS beamformed image in our case), over the normalized autocorrelation 
function (higher than 3 dB) of the images formed by using the other aforementioned BE methods (MV, 
multi-beam Capon, lAA, BP, and LS). Note that a value of RG > 1 needs to be achieved for achieving 
a better resolution than DAS beamformer. 


V. Results and Discussion 

A. Individual point reflectors 

With this simulation we evaluate the potential of the proposed methods in sparse mediums. The resulted 
beamformed images are illustrated in the Fig. The result of DAS BF is shown in Fig. |^a). Using the 
MV BF, the lateral resolution is improved compared with DAS, lAA, and LS BF (see Fig. |^b)), and it is 
comparable with the result of multi-beam Capon BF, Fig. [^c). Concerning the lAA beamformed result, 
as stated in [ fT^ , it gives better point-target resolvability than DAS, Fig. |^d). The proposed BP BF have 
the best resolution of the point-like reflectors, being able to perfectly detect the 5 reflectors, by obtaining 
the most narrower mainlobes, due to the fact that BP results in a sparse representation of the beamformed 
signals. Fig. ^e). As expected, LS beamformer results in solutions that tend to be smooth and regular, as 
in the Fig. |3;f). 

Fig. 1^ presents the lateral profiles of the compared BF methods at 65 mm. We can observe that multi¬ 
beam Capon and MV are comparable in terms of lateral profiles, but MV offers better delimitation of the 
two points. As observed, BP BF outperforms the other BF methods, being able to perfectly resolve the 
two points, suppressing also the sidelobes. Finally, LS BF gives the smoothest result. 
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Fig. 3: (a) DAS, (b) MV, (c) multi-beam Capon, (d) lAA, (e) BP, and (f) LS BF results of the simulation 
of individual point scatterers. 


Lateral profiles at depth 65 mm 



-4.5 -3 -1.5 0 1.5 3 4.5 

Lateral distance [mm] 


Fig. 4: Lateral profiles at 65 mm depth of the point reflectors represented in Fig. 


B. Simulated hypoechoic cyst 

The BF results of a hypoechoic cyst in a speckle pattern are shown in the Fig. We have highlighted 
with white circle the true borders of the cyst, in order to show the accuracy of the proposed methods 
regarding the dimensionality of the scanned structures. 

The image quality metrics are detailed in Table |I^ To calculate the CNR, we have considered the region 
inside the hypoechoic cyst (the black region), and the region Ri inside the homogeneous speckle, at 
the same depth and with same dimension as the region R 2 , as suggested in Q. The SNR was computed 
for Ri. For calculating RG, the whole image was considered. As expected, with DAS the cyst appears 
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Fig. 5: (a) DAS, (b) MV, (c) multi-beam Capon, (d) lAA, (e) BP, and (f) LS BF results of the hypoechoic 
cyst simulation. 


TABLE II: CNR, SNR, and RG values for the simulated phantom in Fig. 


BF Method 

CNR 

SNR 

RG 

DAS 

4.8 

0.4 

1 

MV 

5.3 

0.61 

3.64 

multi-beam Capon 

5.4 

0.58 

4.87 

lAA 

3.5 

0.63 

3.57 

BP 

6.5 

0.62 

8.72 

LS 

7.4 

0.68 

2.65 
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Fig. 6: Lateral profiles at 80 mm depth of the cyst phantom represented in Fig. 


more narrow due to the low resolution and its low capability of resolving cyst-like structures inside the 
speckle pattern, Fig. |^a). By using MV, we slightly increase the contrast and the resolution in the final 
image compared with DAS, the dimension of the cyst being closer to its real dimension, as shown in 
Fig. I^b). Better resolution is obtained when the multi-beam Capon is used, the RG being increased by a 
factor of almost 1.4. The improve in resolution can be observed also in the delimitation of the cyst region 
compared with the white circle that represents the real dimension of the cyst, see Fig. |^c). Compared 
with DAS, lAA increases the resolution of the beamformed image, but not as much as MV or multi-beam 
Capon, Fig. |^d). However, a contrast degradation can be observed from Table Finally, the proposed 
methods are reflecting more correctly the real dimension of the cyst, especially when using the BP BF, 
Fig. I^e), this being in concordance with the high increase in resolution (with a factor of 2 compared with 
MV) and contrast. As expected, LS tends to favor continuity and smoothness, especially when dealing 
with the speckle pattern (see Fig. |^f)), the gain in resolution being less important. However, even so, it 
is more precise in reflecting the dimensionality of the cyst. Note that in terms of SNR, in comparison 
with DAS, all the other beamformers give better SNR, the best improvement being obtained with LS BF 
which is also outperforming the other beamformers in terms of contrast. 

Fig. [^presents the lateral profiles of the results presented in the Fig.|^ where the previous observations 
are confirmed. The curves in Fig. |^are computed by averaging 15 lateral profiles around depth 80 mm. 
The proposed methods, BP and LS, have larger mainlobes than the other BF methods, that correspond to 
the true dimension of the hypoechoic cyst. We can also confirm the increase in contrast presented in the 
Table [H] in case of BP and LS BF approaches. 

Fig .^presents the variation of CNR and SNR parameters function of A hyperparameter. We can observe 
that a favorable compromise between CNR and SNR is reached when A = 0.5. The value of CNR can be 
improved by increasing the value of A. For example, when A = 0.9, CNR= 6.61, but the value of SNR 
is reduced, SNR= 0.55. Similarly, decreasing the value of A will increase the value of SNR, while losing 
in CNR. 


C. Simulated cardiac image 

The results of beamforming the cardiac medium are shown in the Fig. With this example, we are 
interested in visualizing the LV region (hypoechoic), which is surrounded by the hyperechoic regions 
containing the anterior and posterior walls of the heart as well as the septum. The small echoic regions 
inside the LV region are the papillary muscles, that due to the low contrast and resolution of the DAS 
beamformed image are hard to be distinguished. Fig. [^a). A better visualization of the walls is obtained 
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Fig. 7: The variation of (a) CNR and (b) SNR versus A when BP method was applied to the hypoechoic 
cyst simulation. 



E 

E 



Fig. 8: (a) DAS, (b) MV, (c) multi-beam Capon, (d) lAA, (e) BP, and (f) LS BF results of the ultrarealistic 
simulation of a cardiac image. 
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TABLE III: CNR, SNR, and RG values for the simulated US cardiac beamformed images in Fig. 


BF Method 

CNR 

SNR 

RG 

DAS 

1.12 

0.47 

1 

MV 

0.90 

0.56 

3.56 

multi-beam Capon 

0.61 

0.54 

4.23 

lAA 

1.30 

0.62 

5.3 

BP 

1.45 

1.75 

9.89 

LS 

1.55 

1.88 

2.06 


TABLE IV: CNR, SNR, RG, and computational time values for the experimental carotid beamformed 
images from Fig. 


BF Method 

CNR 

SNR 

RG 

Computation time [s] 

DAS 

1.84 

1.46 

1 1 

4.5552 

MV 

2.24 

1.43 

1.25 1 

122 

multi-beam Capon 

1.32 

1.49 

1.25 1 

368 

lAA 

1.48 

1.47 

1.34 1 

8.9266 

BP 

1.85 

1.49 

1.48 1 

60.4320 

LS 

1.94 

1.55 

1.17 1 

8.8692 


with MV (Fig. [^b)) and multi-beam Capon (Fig. [^c)), resulting also in an improved resolution, confirmed 
with a higher RG value (see Table ra For the calculation of the CNR, we considered Ri the region 
inside the white square, situated at approximately 18 mm (laterally) and around 55 mm (axially), while 
is delimited by the black square, around -20 mm (laterally) and 55 mm (axially). To compute SNR, 
the Ri region was considered. 

An interesting observation is that the value of the CNR in the case of MV and multi-beam Capon 
is not improved compared with DAS. This is explained by the results in Q, where it has been shown 
that the improvement of the contrast directly depends on the high definition of the regions (the LV, the 
septum, and the walls in our example). Since the amplitude of the reflectors from the walls and septum 
are not so high compared with the region of LV that contains speckle, the contrast of the final image 
is affected. However, lAA improve both the contrast and the resolution of the image, presenting more 
defined regions, as shown in Fig. [^d). Yet, the best improvement of the resolution is obtained when we 
promote Laplacian BF solutions, with BP BF, see Fig. [^e), resulting in an improvement by a factor of 2 
in RG compared with MV and multi-beam Capon, and by a factor of almost 10 compared with DAS. Of 
course, as expected, LS BF is highly improving the contrast and the SNR of the resulted image, while 
the RG is lower than when using the other BF approaches. Fig. [^f). 


D. In vivo data: carotid 

Fig. 1^ presents the BF results of the studied beamformers, and in the Table we calculated their 
corresponding CNR, SNR, RG, and computational time values. In this example, the carotid is placed 
between 8 and 15 mm in the axial direction. In this region, the interior of the carotid artery is the hypo- 
echoic structure surrounded by the arterial walls (which are hyper-echoic). To calculate the CNR, we have 
considered region R 2 inside the carotide (the white rectangle positioned at 0 mm laterally), and the region 
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Fig. 9: (a) DAS, (b) MV, (c) multi-beam Capon, (d) lAA, (e) BP, and (f) LS BF results of experimental 
carotid data. 


Ri inside the region of speckle (the black rectangle positioned at 0 mm laterally). The SNR for Ri was 
computed. 

As observed, by using DAS BF is hard to distinguish between the interior of the carotid and its walls. 
Fig. I^a). This can be also explained by the fact that DAS BF result represents the lower RG. A better 
visualization of the structures of interest are obtained with MV and multi-beam Capon, that have similar 
RG values. However, the contrast of the MV beamformed image is better, increasing the value of CNR 
by a factor of 1 compared to multi-beam Capon. We can observe that multi-beam Capon is clearly 
defining the region inside the carotid, by reducing the level of speckle inside it, cf. Fig. |^c). The lAA 
beamformed image is comparable with the one of multi-beam Capon, but it conserves better the speckle 
inside the carotid, offering a better resolution and a better contrast of the image. With the proposed 
approaches, however we are able to better distinguish the interior of the carotid artery, as well as its 
walls, with a high gain in contrast and resolution resulted by applying BP BF. A loss in resolution can be 
observed when using LS, compared with BP, due to the use of the f 2 -norm regularization. Note that, due 
to the formulation of the proposed direct model @ that includes an additive noise, the proposed method 
is intrinsically denoising the signal (e.g., the noise inside the carotid is reduced) through the inversion 
process (see Table [TVl). The denoising effect obtained by our BF approach does not suffer from any spatial 
resolution loss, as it could be the case if the raw data or the beamformed images were low-pass filtered. 

Regarding the computational time, note that it is highly dependent on the length of the acquired raw 
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Fig. 10: (a) DAS, (b) BP, and (c) LS BF results of in vivo helthy thyroid data. 


TABLE V: CNR, SNR, and RG values for the in vivo healthy thyroidal beamformed images from Fig. 


10 


BF Method 

CNR 

SNR 

RG 

DAS 

0.55 

0.22 

1 

BP 

1.13 

0.32 

3 

LS 

1.56 

0.64 

2.5 


data. For this case, the number of ranges was around N = 2500, and the proposed approaches were applied 
without a previous decimation of the raw data. Of course, the standard parallel computing methods could 
additionally improve the computational complexity, since the BF process is done for each lateral scanline. 
All the discussed methods were implemented with Matlab R2013b, on an Intel iV 2600 CPU working at 
3.40GHz. Note that even so, LS BF is approaching the time capabilities of DAS, being just twice slower 
than DAS. Moreover, BP is also faster than MV. Thus, by using the discussed techniques for improving 
the computational expense, makes the two proposed methods good candidates for real-time applications. 


E. In vivo data: healthy thyroid 

Fig. [represents the beamforming results of healthy thyroid data. The thyroid (echoic region) is situated 


between the trachea and the carotid artery (laterally, between -20 mm and 30 mm approximately). Fig. 10 


(a) illustrates the result obtained with DAS BF. As expected, the contrast of the image is low and it is 
hard to distinguish the thyroid structure from the trachea, especially in the upper-left part of the thyroid. 
However, when BP (Fig. [T^ (b)) and LS Fig. [T^ (c) are used, the thyroid region is easy to be identified, 
and the contrast of the image is increased. 

The values of CNR, SNR, and RG are depicted in the Table jVj To compute CNR we considered region 
inside the thyroid (the black circle positioned at approximately —10 mm laterally), and the region 
R\ inside trachea (the white circle positioned at approximately —40 mm laterally). The SNR for R\ was 
computed. We can observe that the best values of the CNR and SNR are obtained when LS method was 






















17 



(c) 

Fig. 11: (a) DAS, (b) BP, and (c) LS BF results of in vivo thyroid data with tumor. 

TABLE VI: CNR, SNR, and RG values for the in vivo thyroidal beamformed images from Fig. [IT] 


BF Method 

CNR 

SNR 

RG 

DAS 

0.71 

0.62 

1 

BP 

1.16 

0.79 

2.9 

LS 

1.32 

0.86 

1.5 


applied, the thyroid region being obvious to be discerned. The boundaries of the carotid artery are also 


well defined. Fig. 10 (c). 


F. In vivo data: thyroid with tumor 


The beamformed results of the thyroid data with tumor are presented in the Fig. 11 The malignant 


tumor with an irregular structure can be seen between the left lobe of the thyroid (the hyper-echoic 
structure situated near the trachea) and the carotid artery (the hypo-echoic circular structure with the 
center at approximately 33 mm (axially) and 40 mm (laterally). We can observe that, contrarily to DAS 
beamformed image, where the tumor is hard to be distinguished (see Fig. [TT] (a)), both our methods 
improve the visualization of the main structures, enhancing the edges of the tumor. The values of CNR, 


SNR, and RG are depicted in the Table where a gain in resolution with a factor of almost 3 can be 
observed when using BP, compared with DAS, while with LS we obtain a higher improvement in contrast 
and SNR than with BP BF. CNR was computed by considering region R 2 inside the tumor (the black 
circle positioned at 0 mm laterally), and the region inside the left lobe of the thyroid (the white circle 
positioned at approximately —20 mm laterally). The SNR for Ri was computed. 


VI. Conclusion 

We have presented a new BF approach in US medical imaging, that solves a regularized inverse problem 
based on a linear model relating, for each depth, the US reflected data to the signal of interest. Contrarily 
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to existing techniques that use adaptive or non-adaptive weights to average the raw data in order to 
form RF lines, we directly recover, for each depth, the desired signals using Laplacian or Gaussian 
statistical assumptions. The proposed regularization-based BF allows us to take advantage of the beamspace 
processing that enables to highly reduce the number of US transmissions (by a factor of five in our 
examples), while improving the quality of the beamformed images compared to four existing beamformers. 
Multiple simulated and experimental examples were presented, that compare our approach with DAS, MV, 
multi-beam Capon, and lAA beamformers. We showed that our BF approaches, based on Laplacian and 
Gaussian prior information, although based on the same model, are complementary in terms of result 
quality. Thus, Laplacian statistics are favoring sparse results while the Gaussian law is offering more 
regular and smooth images. We also proved through resolution gain, CNR and SNR image quality metrics, 
that we obtained an important gain in spatial resolution and/or in contrast, while maintaining a reasonable 
computational time compared to other existing techniques. As future work, we will consider other statistical 
assumptions, such as the generalized Gaussian distribution, resulting in £p-norm minimization with the 
parameter p between 0 and 2. Following the choice of p, this should guarantee a better compromise 
between the gain in contrast and the improve of the spatial resolution (e.g. p6| , | |47] |, and Another 
way to obtain this compromise could be to combine the two regularization terms (through Laplacian and 
Gaussian statistics) used in our approach, resulting in an elastic net regularization (see e.g., [ |25] | and [|48]|). 

Another interesting perspective offered by our BF direct linear model is the possibility to combine it with 
existing post-processing techniques, aiming to enhance the quality of US images, such as deconvolution 
or super-resolution. 


Appendix A 

Minimum variance beameorming 

MV (or Capon filter) BF |j^ consists in minimizing the array output power by maintaining a unit gain 
at the focal point. It adaptively calculates the weights, by solving: 


mmw^RkW, such that = 1, (19) 

W 


with the analytical solution: 


Wmv 




( 20 ) 


where = E[y^y^] is the covariance matrix of y^ and 1 is a length M column-vector of ones. These 
weights are then used to calculate the desired RF beamformed lines in a similar way as with DAS. In 
practice, R^ is unavailable and the estimated covariance matrix R^ is used as alternative, derived from 
L received samples: 


L 

l=\ 

Since the received focused raw data is coherent, several methods were proposed to decorrelate the data 
as much as possible: subaperture (or subarray) averaging (also called spatial smoothing), time averaging, 
and diagonal loading significantly improve the standard MV BF (e.g., [|50|). 

Appendix B 

Beamspace beamforming 

Starting from the MV BF method presented in Section named element-space based Capon (ES- 
Capon) in pOj , Nilsen et. al. proposed a beamspace beamformer (BS-Capon) that allowed reducing the 
computational complexity of the MV BF by a ratio of 3. Basically, they reduce the size of the covariance 
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matrix in ( |^ by replacing it with a smaller covariance matrix of orthogonal beams. The expression of 
the orthogonal beams is detailed in p0||, and the beamspace transformation is expressed as follows: 


VkBS = BVk. 

where B = \bi, - ■ ■ , buY is the M x M Butler matrix whose elements are defined as: 


bmn rr-;: 

VM 




( 22 ) 


(23) 


B is an unitary matrix (BB^ = B^B = Imxm), equivalent to an M-point discrete Fourier transform 
(DFT) matrix. ImxM is the identity matrix of size M x M. 

The transformation in ( |22l ) is applied to all signals and weights vectors in the element-space (ES) to 
find their beamspaced version. Therefore, the weights of ES-Capon BF are formed by solving: 


tH: 


min w^gRss'iVBS, such that 


WBS 


w^s^i = 1 - 


The solution of (24) is: 


wbs = 


YR-B\e^ ’ 


(24) 


(25) 


where Rbs = is the covariance matrix of y^gg is a M x 1 vector having the value 

1 in the m-th position and zero in all other positions. Finally, we can state that BS-Capon BF can be seen 


as the description of the Capon filter from ( |20| ) in the Fourier domain. 

As stated before in this section, in the case of DAS, MV, and BS-Capon BF, the final RF US image is 
a collection of RF beamformed lines, each of which being the result of beamforming the raw RF signals 
coming from a focused wave emission in the direction 9^, k ^ {1,... ,K}, using the M elements of the 
transducer. 

The most commonly used visualization mode in ultrasound medical imaging is the B (brightness)-mode. 
It is obtained by applying envelope detection and log-compression techniques to each beamformed RF 
line. Finally all the RF lines are juxtaposed in the lateral direction to form the final 2D US image, as 
shown in the Fig. 


Appendix C 

Multi-beam Capon beameorming 

Jensen et. al. used beamspace processing for reducing the dimensionality of the data, and proposed a 
new approach of Capon BF, called multi-beam Capon BF ||TT|. For more convenience, let us briefly recall 
their approach. 

For a given range n, let us select its corresponding lateral scanline, as illustrated in Fig. Since the 
signals ykgg have been focused in axial direction (by applying time delays) before being beamformed, we 
just need to compensate the phase-shifts based on the distances from the samples of the lateral scanline 
(equivalent, in our case, with K, the number of beam directions). 

The compensation of the phase-shifts, with M = 0, • • • , M — 1 is depicted in Fig. [T^ assuming 

that the time-compensated data reaches the elements at angle Ok- We consider the first elements as 
reference, so its phase-shift is 0. Thus, based on the far-field assumption, we can formulate the complex 
exponential version of the manifold vector for a given direction k, that corresponds to the incident angle 
6k, as (see e.g. dD and [ |51| |): 


a0^ = [1 . . . g-i(M-l)7rsm(0fe)jr^ 


(26) 
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Fig. 12: Phase shift compensation of the focused raw data. 


Thus, by using phase shifts, for focusing along a lateral scanline, contrarily to the matrix Rbs used in 
BS-Capon, for a given range n, the covariance matrix R[n\ will cover all the directions 6k,k = 1, • ■ • ,K. 
Therefore, the weights corresponding to a given direction 9 and a range n will be formed by solving: 


mmw^R[vi\w, such that w^ae^n = 1, 




(27) 


having the solution: 




R [n\a0^r. 


(28) 


'^[n]a0^n 

These weights are then applied to calculate the signal corresponding to a lateral scanline, at a range n. 


Appendix D 

Iterative adaptive approach beamforming 

Based on the beamspace processing technique and on the calculation of the multibeam covariance 
matrix discussed in the Section]^ Jensen et al. applied lAA [ fT^ to US medical imaging, [ fT2| |. Following 
this recent work, a covariance matrix, R[n] based on K potential reflectors placed across a considered 
lateral scanline, was defined as: 


K 

= AbsPAIs, (29) 

k=l 

with VBsb^] ^ the beamspaced time-delayed raw data at a given range n, before applying the 

phase-shift transform. A is the matrix containing the manifold column-vectors defined in ( |26l ), and P a 
diagonal matrix with the elements of along its diagonal. The values of P are then iteratively 

updated and calculated by taking into account the weights corresponding to a lateral scanline, by following 
( |28] ). Finally, P is used to estimate the amplitude of each reflector of the lAA BF result. 

Contrarily of DAS, MV and BS-Capon BF where, to form the final beamformed image, the RF lines are 
juxtaposed in the lateral direction, multi-beam Capon and lAA BF are axially juxtaposing the beamformed 
lateral scanlines to form the final beamformed image. 
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